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Abstract 

The geometry of a light wavefront evolving in the 3-space associated with a post-Newtonian 
relativistic spacetime from a flat wavefront is studied numerically by means of the ray tracing 
method. For a discretization of the bidimensional wavefront the surface fitting technique is used 
to determine the curvature of this surface at each vertex of the mesh. The relationship between 
the curvature of a wavefront and the change of the arrival time at different points on the Earth is 
also numerically discussed. 



I. INTRODUCTION 



The description of the propagation of hght in a gravitational field is even today a central 
problem in the general theory of relativity. The deflection of light rays and time delays of 
electromagnetic signals due to the presence of a gravitational field are phenomena detectable 
with current experimental techniques which allow design new tests for general relativity. In 
this line, Samuel p^j recently proposed a method for the direct measurement of the curva- 
ture of a light wavefront initially fiat, curved when light crosses regions where the gravita- 
tional field is non vanishing. He found a relationship between the differences of arrival time 
recorded at four points on the Earth, measured by employing techniques of very long base 
interferometry and the volume of a parallelepided determined by four points in the curved 
wavefront surface. This surface is described by means of a polynomial approximation of the 
eikonal in a Schwarzschild gravitational field. For more complex gravitational models, such 
as those considered by Klioner and Peip [2], de Felice et al. [3] or Kopeikin and Schafer |3] 
in studies of light propagation in the solar system, the use of numerical methods for the 
determination of the geometry of the wavefront surface would also be required. An ana- 
lytical approach to the relativistic modeling of light propagation has also been developed 
recently by Le Poncin-Lafitte et al. and Teyssandier and Le Poncin-Lafitte |6], where 
they present methods based on Synge's world function and the perturbative series of powers 
of the Newtonian gravitational constant, to determine the post-Minkowskian expansions of 
the time transfer functions. 

Nowadays there are numerous techniques in computational differential geometry which 
allow to analyze geometric properties of surfaces embedded in the ordinary Euclidean space. 
Techniques of this type are widely applied in different areas such as Computational Ge- 
ometry, Computer Vision or Seismology. In one of these methods, developed in works by 
Garimella and Swartz [7] and Cazals and Pouget [S] , the estimation of differential quantities 
is established using a fitting of the local representation of the surface by means of a height 
function given by a Taylor polynomial. A survey of methods for the extraction of quadric 
surfaces from triangular meshes is found in Petitjean [9]. 

In this work we consider a discretization of the wavefront surface, replacing this surface 
by a polyhedral whose faces are equilateral triangles. At initial time, the surface is assumed 
to be fiat and far enough from a gravitational source (say, the Sun) and moving towards 
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this source. We study the deformation of the instantaneous polyhedral representing the 
wavefront when crossing a region in the relativistic 3-space near the Sun due to the bending 
of light rays by the gravitational field. In this study, we apply the ray tracing method with 
initial values on the vertices of the triangular mesh to obtain the corresponding discrete 
surface at each instant of time. Then we apply the techniques given in |7j and [8j to describe 
the wavefront as a surface embedded in the Riemannian 3-space of the post-Newtonian 
formalism of general relativity. For each vertex in the instantaneous mesh we obtain a 
quadric which represents locally the surface by applying the least-squares method to the 
immediate neighboring vertex around the considered point which is represented in normal 
coordinates adapted to the light rays. 

The structure of the paper is as follows: In Section 2 we briefly introduce the basic model 
for the wavefront propagation in the post-Newtonian formalism. In Section 3, we establish 
a discretized model of the wavefront surface by means of a regular triangulation and we 
describe the method employed in this work for the study of the curvature of this surface. 
In Section 4, a numerical estimation of the curvature of the surface is derived using the 
ray tracing method. For the numerical integration of the light ray equation, we use the 
Taylor algorithm implemented by Jorba and Zou [10] which is based on the Taylor series 
method for the integration of ordinary differential equations and which allow the use of high 
order numerical integrators and arbitrary arithmetic accuracy, as is required to describe the 
influence of weak gravitational perturbations on the bending of light rays. Finally, in Section 
5, we apply the method discussed above to study the effect of the wavefront curvature on 
the variation of the arrival time of the light at points on the Earth surface, following the 
model proposed in Samuel's test [Ij. The paper concludes by giving another approach to the 
estimation of the curvature of the wavefront, derived from an approximation of the Wald 
curvature [12] associated with a quadruple of points in the wavefront. 

II. LIGHT PROPAGATION IN A GRAVITATIONAL FIELD 

Let us consider a spacetime (^, g) corresponding to a weak gravitational field and choose 
a coordinate system {(z, ct)} such that the coordinate representation of the metric tensor is 

9ap = Vap + Kp, with Tjan = diag (1, 1, 1, -1). (1) 
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where the coordinate components of the metric deviation hap are given by: 

hah = 2c^^K||z||^"^5aft, hai = —4:C^^ K,\\z\\^^ Za, /;.44 = 2c^^/t || 2 1| (2) 

(Greek indices run from 1 to 4 and Latin indices from 1 to 3) where k := GM represents the 
gravitational constant of a monopolar distribution of matter (say the Sun) located at 2'"(t) 
and c represents the light speed. 

In the post-Newtonian framework one may consider a simultaneity space at each 
coordinate time t. From the fundamental equation of the geometrical optics for the phase 
ipizjt) of an electromagnetic wave jH]: 

.«/3^^_0 (3) 

and Cauchy data 'ip{z, 0) = const, given on a spacelike surface '■= {(^5 0) | (p^z) = 0} C Sq 
one obtains the characteristic hypersurface (light cone) Vl := {{z,t) \ ilj{z,t) = const} C 
^ of the light propagation. The intersection of the characteristic hypersurface and the 
corresponding simultaneity space is the spacelike wavefront at a time t which will be denoted 

by yt-=^r] St. 

An alternative formulation of the problem of light propagation in the spacetime (^, g) 
is based on the determination of the bicharacteristics generated by the isotropic vectors 
k := grad'?/'. From both the equation for the null geodesies, z{t) = [z(t),t) expressed in 
terms of the coordinate time and the isotropy condition g{z, z) = 0, one obtains in the 
post-Newtonian approach to general relativity that, neglecting terms of order 0(c~^), the 
null geodesies of (^, g) must satisfy the equations 

Z'' =y^aiz,Z,t) (4) 

=gapz"z''. (5) 
In ^ the components of the acceleration (fa{z,z,t) are given by (see [13]) 

iPa{z, Z, t) = lc^hu,a " [i^44,t<5fc + hak,t + c{hia,k " hik,a)]z'' 

— (/l44,fe'5r + hak,l — \hkl,a)z^Z^ 

-{c-'Kk,, - \c~%k,t)z^zH\ (6) 



where the first and third terms in the right hand side of (|6]) are of order 0(1), while the 
remaining terms are 0{c~^). 



For initial values (20,20), with 20 G ,5^^ and (20/c, 1), satisfying the condition (|5]), the 
integration of the initial value problem corresponding to Q on an interval [0, T] allows to 
determine the spacelike wavefront S^t- This surface is embedded in the Riemannian 3-space 
(Sy, 7) where the components of the metric tensor are given by 

lah ■■= gab , (7) 

544 

and at the time T the tangent vectors z{T) to the integral curves z{t) are 7-orthogonal to 



III. NUMERICAL DESCRIPTION OF A SPACELIKE BIDIMENSIONAL WAVE- 
FRONT 

Hereafter, we consider the simplest gravitational model generated by a static point mass. 
Let S denote the quotient space of ^ by the global timelike vector field dt associated with 
the global coordinate system used in the post-Newtonian formalism. We will consider a 
region of the wavefront in S described by a coordinate chart {2;}. Further, we assume that 
the Riemannian manifold (<f , 7) is almost fiat and that the metric corresponding to 7 is 
quasi-Cartesian in the chosen coordinates. 



A. Discretization of the initial wavefront 



Given an asymptotically Cartesian coordinate system, we consider a set formed by 
points with coordinates (2:1, ^2, — C) where > is a number large enough so that o5^o niay 
be considered as a fiat surface. The direction determined by the point O* : (0, 0, —Q in ,5^q 
and the center of the Sun O : (0, 0, 0) is perpendicular to 5^^. We will study the geometry of 
a region ^ C .5^^ determined by points P whose Euclidean distances d{P,0*) to the point 
O* satisfy Rq ^ d{P,0*) ^ 2Rq, where Rq is the radius of the Sun. 

For the discretization of the problem we consider, in the first place, the set of points 
(ai, 02, as) G where = {0, Ni, Ni + 1, . . . , N2} and Ni < N2 are two natural numbers 
(see Figure [ifa)). In the point (0,0,0) is excluded. Next, we construct in the complex 
plane a regular triangular mesh whose vertices are located between two hexagons as shown 
in figure mb), and the edges have length i. The inner hexagon has sides of length iNi and 
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(3,0,4) 



(a) 



(b) 




75 76 77 78 79 80 



FIG. 1: Triangulation and enumeration of the initial wavefront. (a) Triplets (ai, 02, 03) G used 
to label the vertex of a hexagonal mesh "V . (b) Regular triangulation of a hexagonal annular 
region in the wavefront. Points (3, 4, 0), (0, 3, 4) and (3,4,0) in Y correspond with points 36,44 
and 50 in ^* respectively. 



the outer hexagon sides of length ■ In this triangulation each vertex is represented by a 
complex number of the set (see jl4j ) 

:= = ai + 021^ + «3^^^ I «!, 02, 0.3 ^ ^ ■, ^ '■= exp(27ri/3)}, (8) 

where i := The vertices (01,02,03) with some of their components equal to A^^i (resp. 

A'"2) are located on the inner (resp. outer) boundary of the mesh Y. We establish an 
enumeration of the vertices as shown in Figure [l|^b), in such a way that the inner vertices 
Zj have subscripts j = 1, . . . , J. Finally, we apply the change of scale: 

— ZT 

f^r\ z^—, (9) 

so that the inner boundary is a hexagon of radius r. In consequence the length of each edge 
in this triangulation is equal to r /Ni. 

The complex plane and the plane may be identified by means of the mapping l : z ^ 
(3fJ(z), '^{z)^ ~C)- Thus one obtains a discretization of the initial region of which we are 
considering here, triangulated with vertices given by L{zj) whose corresponding mesh on 
will also be denoted by Y* . 
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B. Normal coordinates around a point 



Now we assume that at each vertex in "V* a photon with velocity Zq := (0, 0, c) is located. 
The null geodesies equation Q may be written as a first order differential system ii = F{u, t) 
in phase space u := {z, z) which determines a flow in 

z{t) =$t{Zo,Zo), (10) 

in terms of the initial values Zq := z(0) G Zq := z(0). Then, for each time t there 



is a surface S^t image of S^q under the flow (10). At a point z ^ S^i the normal vector 



t(2) coincides with the tangent vector to the curve zit) at that point. Furthermore, before 
reaching the focal points of the beam of light, the triangulation "V* induces a triangulation 
"V on the wavefront S^t whose vertices we enumerate using the same labels used for the 
corresponding vertices in 

In order to simplify the description of the geometry of the surface on a neighborhood 
of a point P G we use a 'y-orthonormal reference frame {e.^^^^ centered on that point, 
where one of its vectors, say 63, is parallel to the vector n(P) := t(P)/ a/7(t, t) tangent 
to the ray passing through that point. Let {y-^ jl^^ be a normal coordinate system with pole 
at the point P and associated normal reference frame {ej}?^^ . 

By using the classic formulae of Riemannian geometry (see [15] §18, and p5|), the coor- 
dinate transformation ^ z*, from normal to post-Newtonian coordinates, is determined 
by 

v' = (A-^)U^" - z^o + Kr^e)o(^' - - 4)), (11) 

where A is a non-singular constant matrix and (r^cjo ^td^) Christoffel symbols 

at the point P. By neglecting terms of order higher than F^^, the inverse transformation 
may be approximated by 

z' = < + M{y'-K^k)oy'y')^ (12) 

whose corresponding Jacobian determinant is 

^=K(s}-m^'). (13) 

In normal coordinates a metric tensor 7 on the space S' is determined from 7 by 

= a^^"- ^^^^ 
7 



and, as it is well known, at point P the tensor is reduced to Sij and the associated 
Christoffel symbols at this point are F*^ = 0. 



C. Local approximation of the wavefront 

From the triangulation 'f^* of the initial wavefront the ray tracing method furnishes a 
discrete surface determined by the mesh Y. To compute differential magnitudes of the 
wavefront surface corresponding to this mesh one needs to define a discrete neighborhood of 
each vertex in For this, we consider firstly at each inner vertex z* G = 1, . . . , J a 

neighborhood (named hereafter 1-ring in the terminology of computational geometry, e.g. 
[18] ) formed by the six vertices z*^ closest to z*: 

[2;*; 4Jj=i,...,j,fc=o,...,5, where z*^ := Zj + exp{kni/3), (15) 



Then, for each 1-ring in (15) one may determine on the mesh Y a corresponding 1-ring 



formed by the image of the points z*, z*^ G ^0 under the flow (10) 



[z,;z,^]:=[<P,iz*);<P,{zl)]. (16) 

In a neighborhood of a point Zj the wavefront can be approximated by a height function 
on the orthogonal plane to 63 determined by means of a least-squares fitting of the data 



(16) expressed in normal coordinates {y^}^^^. As a model for this surface we chose a quadric 
passing through the coordinate origin whose gradient at this point is parallel to 63, 

y' = fiy\y') ■= k^iiy'f + a2yV + la.iy'f , (17) 

where ai,a2 and are the indeterminate coefficients to be obtained by the least squares 
method. 



The quadric (17) provides a surface S^j that approximates the surface S^t in a neigh- 
borhood of the point zj and is defined in parametric form y^ = y'^{x^), with A = 1,2, 
as 

y^=x\ y^ = x\ y^ = f{x\x^). (18) 

In coordinates (x^,a;^), the metric 7 on ^ induces a metric on whose associated metric 
tensor g has the form 

dy^ dy^ 
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Moreover, on the tangent plane to S^t at Zj one defines a tensor B associated with the normal 
n at that point as: 

B : T,^.^ X T^.S^ - (T.^,^)^ 5(9^, ds) = g^ri (20) 

Therefore, if {vi, ^2} represents an orthonormal basis for the vector space T^.y consisting 
of eigenvectors of B with associated eigenvalues Ai and A2, and II denotes the second 
fundamental form, then the difference of sectional curvatures K and K associated with 
the plane generated by {vi,i;2}, in and S respectively, is given by a generahzed Gauss 
formula (see [19], p. 131) 

ir,ei := K{vi, V2) - Kivu V2) = A1A2 (21) 

which is named (see [15]) relative sectional curvature -ft'rei, whereas the mean curvature is 
determined by half the trace of 11: 

if=l(Ai + A2) (22) 

IV. NUMERICAL ESTIMATION OF THE CURVATURE OF A WAVEFRONT 
A. Integration of the equations of light rays 

Here we deal with the problem of the numerical integration of the initial value problem 
corresponding to (|4])-(j5]). The equation Q for the light propagation in the gravitational field 
generated by a static material point, may be rewritten in terms of new variables defined as: 

u := {ui, . . . ,Ue), with Ui := Zi and := Zi, (2 = 1,2,3), (23) 

obtaining the following first order differential system: 

ill = 

U2 = U5, 



% = Me, 



U4 



K C^Ui — ?)Uiu\ — 4:U2U4U5 — 4:U3U4Uq + Uiuf + Uiu\ ^24) 



C2 {ul + Ul+ m2)3/2 

K C^Z2 — 4M1U4U5 — ?>U2U\ — Au^U^Uq + M2M4 + ^2^6 

[uj + Ul + Mi)3/2 
K C^Z3 — AUiU^Uq — AU2U^Uq — 3U3Ug + + 

{ul + Ul + Mi)3/2 • 
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The solution u(t) must satisfy the constraint ([s]) which, in the notation (23), may be ex- 
pressed as 

2k 

F{u) := ul + ul + ul-c^+ .{ul + ul + ul + c') = 0. (25) 

c^y + ^2 + ^3 



To obtain a numerical solution of the initial value problem given by (24) and the initial 
data u{0) = (2(0), 2(0)) we use the Taylor integrator developed by Jorba and Zhou [TO] , 
based on the classic Taylor series method for ordinary differential equations. In this method 



a Taylor expansion of the vector field ii defined in (24) is made using techniques of automatic 
differentiation to obtain the corresponding Taylor coefficients. The Taylor integrator allows 
the control of both the order and the step size employed in the method. Furthermore, the 
Taylor integrator is implemented so that one may use extended precision arithmetic for the 
highly accurate computation required in this problem. 

From now on we use normalized units taking the radius and the mass of the Sun as units 
of length and mass, respectively. Then, the initial values corresponding to a photon initially 
located at 100 astronomical units from the Sun are given by 

Uo := u{0) = (1.0, 0.0, -21494.6550, 0.0, 0.0, 0.4307). (26) 

For the integration over a period of time of 50400 seconds, using a precision of 120 binary 
digits and a tolerance Tol = 1. E— 20, the Taylor algorithm gives the solution in terms 
of Taylor polynomials of degree twenty-four. As a test of the accuracy of the method, in 
Figure [2] it is shown the behavior of the function F(u) appearing in the isotropy constraint 
(|5|. In this fi gure we see that at the arrival point, the isotropy constraint on the tangent 
vector to the light ray is not satisfied with the same degree of accuracy required at the 
initial position, reaching this deviation its maximum value at the instant when the photon 
is nearest to the Sun. 

In order to obtain from the numerical solution given by Taylor another value Un 



satisfying the condition (25) and such that the function — iinW reaches a minimum, we 
will apply at each step in the Taylor integrator the method of standard projection [20] to 
project Un on the manifold F{u) = 0. This leads to a constrained extremum problem with 
a Lagrangian function L{u, A) := ^\\u — u\\^ — F{un)^X, where A is a Lagrange multiplier. 
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FIG. 2: Semilogarithmic representation of the difference F{un) —F{uq) along the numerical solution 
of the light ray equation obtained by means of the Taylor algorithm. The maximum value of AF 
is reached at the time when the photon is nearest to the Sun. 



The necessary condition of extremum and the constraint condition leads to 

Un = Un + VF(U„)'^A 

= F{un) 



and replacing (27) in (28) the following nonlinear equation for A is obtained 



(27) 
(28) 



(29) 



This equation may be solved by applying the simplified Newton method. The projection 
stage in the algorithm spent a 3% of the 0.06 seconds-CPU employed by an Intel Core 2 Quad 
processor to determine the trajectory of a photon in the time interval we are considering. 



B. Deformation of a wavefront by a static gravitational field 

Now we apply the method described in Section 3 to that region of a wavefront propagating 
along a tubular neighborhood around the Oz'^-axis and whose radius is 2Rq. Initially, the 
wave surface is flat, perpendicular to the Oz^-axis and with position and velocity given 



in (26). The wavefront S^t is then determined after a trip of 101 astronomical units. 



1. Curvature at a point 

Firstly, we obtain an estimation of the curvature at a point in by using a sequence 
of 1-rings with decreasing radii until a stable value of the curvature wavefront at that 
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FIG. 3: 1-ring with vertices [Vq, Vi, . . . ,Vq] and corresponding least-squares quadric fitting of these 
vertices represented in normal coordinates centered at Vq. 

point is reached. Let us consider the 1-ring of radius r and centered at the point V^'^^ := 



(1.0, 0.0, —21494.6550), determined by the three first components in (26) and whose vertices 
are given by 

Vj'^+il := (1 + r cos(A;7r/3), sin(A;7r/3), -21494.6550) , (A; = 0, ... 5). (30) 
Now, we apply the Taylor algorithm for the values of tolerance and arithmetical precision 



pointed out in Subsection IV A , to determine the images l^W of the vertices V^'"\ k = 0, . . , 6, 
by solving (24) and taking for all vertices the same initial velocity (0.0,0.0,0.4307) (in 
the normalized units we are considering). After a time T = 54400 seconds, the 1-ring 
[y[o] - _ ^ y[6]j provides a discretized neighborhood of the point V^^l G (see Figure |3]). 



Applying the scheme developed in Subsection III C we obtain that the mean and the 
relative total curvatures at point 1^'°^ take, for both the different values of the radius r 



and the tolerances in the Taylor algorithm, the values shown in Table IV B 1 where one 
observes that for values Tol = l.E— 20 and r = Rq/50 the three first significant decimal 
numbers are correct. 

Table: Variation of the relative total curvature and mean curvature with respect to both 
the radius of the 1-ring and the tolerance used in the numerical integration of the ray 
equation. 
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,2. Wavefront surface 

To obtain an estimation of the curvature of a region of the wavefront propagating along 
the Oz^-axis, we apply the ray tracing method with initial values on the wavefront surface 



described at the beginning of Subsection IV B We use a regular discretization of the wave 



front surface as that described in Subsection |III A[ Furthermore, taking into account the 



results shown in Table IV B 1 we choose the length of the edges in the corresponding mesh 
equal to i?0/5O. 

In the numerical model we are studying, one assumes that the Sun is a point and we 
consider a hexagonal annular region on the initial wavefront similar to that shown in Fig- 
ure [T|^b) where the inner and outer hexagons have radii of lengths Rq/25 and 2i?0 respec- 
tively. Therefore a number of := 28044 vertices is required. The Taylor algorithm with 
projection in a time interval [0,T] (the time required to run a path of length equal to 101 
AU) is applied to each photon located at an initial vertex; the CPU-time employed to carry 
out this computation is of 1544 seconds. 



For the numerical solution Un, n = 1, . . . , A^, of (24), both the mean and relative-total 



curvatures at each inner vertex of the mesh on the surface are computed by applying 



the method described in Section III, schematized in pseudocode in the next Table. 
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(z* z*^ 



1,...N 



1 1 rays tracing 



// see Eq. (11) 



// see Eq. (16) 



// see Eq. (18) 



Data: : 
for n = 1 . . . do 

Un ■■= Taylor(t, w;) 

:= NormalCoordinates(z,„) 
for i = ... 6 do 

:= Ring(z/„) 

end 

(01,02,03) = LeastSquares(y„J 
lAB{Xn) := Metric(ai, 02,03, x„ 
5 = SecondFundamentalForm(a3„) // see Eq. (20) 
(Ai, A2) = Diagonalize(i?) 
{Kj-ei, H) = Curvature(Ai, A2) 
end 

In Figure |4| the surface at the time when the wavefront arrives at the Earth is shown 
using a gray-scale to represent the mean curvature (note we have used a different scale on 
the Oz^-axis). One sees in this figure that the absolute value of the mean curvature function 
defined on increases as the distance between the photon and the 02;'^-axis decreases. 



// see Eq. (19) 



// see Eqs. (21 



22) 



x10 

6^ 




X 10 



-2 




FIG. 4: Wavefront surface and relative total curvature (gray scale) deformed by a spherical gravi- 
tational field (a different scale is used for the vertical axis). 



14 



C. Curvature of a wavefront in the PPN formalism 

To derive the deflection angle for light rays, instead of assuming the validity of general 
relativity, one may consider a more general expression for the metric generated by a spherical 
central body that is valid for different gravitational metric theories. In the parameterized 
post-Newtonian formalism, the expression of a spherically symmetric metric, written at 
the linearized order we are considering here, contains one parameter 7 which is usually 
interpreted as a modification of the curvature of the space. In the parameterized post- 
Newtonian formalism the total relativistic deflection angle of light rays passing near the 
limb of the Sun is given approximately as (see [21]): 

A(/. ~ 2(1 + 7)1".75. (31) 

Using very long baseline interferometry techniques, one may obtain high precision general 
relativistic measurements for the deflection of radio signals from distant radio sources with 
an accuracy at the 0.02 percent level. In [22j an estimation of 0.9998 ± 0.0004 is given for 
the post-Newtonian parameter 7 . 

Here we consider a gravitational field depending on the Eddington parameter 7 and 
adapt the numerical method described above to determine the dependence of the wavefront 
curvature at a point, corresponding to a light ray grazing the Sun with respect to this 
parameter. For the numerical discussion, we take a gravitational field in which only the 
dominant terms in the post-Newtonian metric are included, so that the metric deviations 
may be written as 

hab = '2c~'^K'-f\\z\\~^5ab, = 0, /l44 = 2c"^k|| Z || (32) 

We take values of 7 in the interval [0.9992, 1.0012] in which the experimental values given in 
[22] are included and consider eleven nodes in the numerical code developed above obtaining 
the results shown in Figure [5] for the curvature of the light wavefront at a point near the 
Earth for a light ray grazing the Sun. In this figure one observes a linear dependence of 
the total and mean curvature with respect to the parameter 7. For these values of 7 the 
variations of i^rei and H are in the second and third significant decimal digit respectively. 
In consequence, the numerical treatment we applied gives account of the variations of the 
geometric model when a change of the PPN Eddington parameter is made. 
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0.9991 0.9996 1 1.0004 1.0008 

Y 



FIG. 5: Relative total curvature i^rei (a) and mean curvature H (b) versus post-Newtonian pa- 
rameter 7 for a wavefront numerically determined. 



V. VARIATION OF THE TIME OF ARRIVAL AND CURVATURE OF THE 
WAVEFRONT 

In this section we consider the problem treated by Samuel [1] where the curvature of a 
wavefront in a gravitational field is detected by measuring the arrival times T^, a = 0, 1, 2, 3, 
at four receiving stations located at points Ea on an Earth hemisphere. The arrival time 
differences between these stations depend on the curvature of the wavefront and are related 

¥ 

to the volume of a parallelepiped formed by the vectors EqEj — c(Tj — To)nj, j = 1,2,3, 
where n := (0,0, 1). Samuel proposed the measurement of this non-zero volume as a new 
test of the general relativity. 



A. Volume of a tetrahedron with vertices on the wavefront surface 



Here we compute the volume and the arrival time differences associated with a simplex 
determined by four points on the wavefront numerically determined above. For this, we 
consider a light ray reaching the reference station Eq (see Figure |6]) at the instant of time 
To and whose tangent vector in this point is Uq. Let Pq be the position of this photon at 
a previous instant Tq such that c{Tq — Tq ) = 101 AU. This position may be obtained by 



a backward integration of the equations of motion (24) with initial values {Pq, — cno). On 
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the plane determined by (Po,nQ), where is the tangent vector to the ray at the point 



Pq, we consider the 1-ring [Pq; P*], j = 1, ... ,6, centered at Pq. Then, from (24) and the 
initial values {P*,cnl) we determine the 1-ring [Qo',Qj] image of [PQ;Pj'] by the flow (10) 
integrating again over the interval [0, Tq — Tq]. 

In the numerical model we are considering, the coordinates of Pq are {2Rq, 0, 1 AU) and 
no = (0, 0, 1). The original position Pq and direction Uq of the photon are then (expressed 
here using only three decimal digits): 

Pq := (1.909, 0.0, -21494.224), n* = (0.183 E-5, 0.0, 0.431). (33) 

We choose as vertices P* the points 

P; = P(pj,e.,^)(r cos(7rj73),r sen(7rj73), -WOAU) J = 1, ... ,6 (34) 

where P(p^*^e2,i?) represents a rotation around axis (Pq , 62) through an angle ^ := e3-nQ/||nQ||. 



Following the same method employed in Subsection IV B 2 we perform a least-squares 
fitting of the data {Qa}a=o to get an approximation of the wavefront surface in a neighbor- 
hood of Qo given by a quadratic surface Q which, in the normal coordinates corresponding 
to ray (Pq, no), takes the form 

Q{z) : Z3 = -0.1060 X IQ-h^ + 0.6 x 10~^^ZiZ2 + 0.1065 x 10-^;z^. (35) 

Now, on the Earth surface we choose four points Ea with geocentric coordinates 
(0,0,-Pe),(^e,0,0),(^Pe4^®'0) ^nd (^Pe,-^^e,0) respectively, which will be 
transformed to normal coordinates. Assuming that the geometry of the 3-space in the 
vicinity of the Earth is Euclidean, one may determine the points Qa ^ ^ whose distances 
to the corresponding stations Ea are minima. 

Once the points Qa & ^ corresponding to the images of the stations Ea on the Earth 
are determined, we may compute the differences Tq = dist((5a, QD/c between the arrival 
times measured at each station under the assumptions that the wavefront surface is either a 
curved surface ^ or a plane 1^ determined by the pair (Pq, no). These differences are given, 
for the data provided above, by 

To = 0.97872481493, n = 0.99999999032, = = 0.99999999021, (36) 
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whereas the distances kj := dist(Qj, Qj), i,j = 0, 1, 2, 3, between the projected points Qa on 
Q are (in the normahzed units we are using) 



^01 - 
/l2 



0.9164245 E-2, 
= 0.4743770 E-2, 



lo2 - 
/l3 



0.9164249 E-2, /qs 
= 0.4743770 E-2 I23 



0.9164249 E-2, 
0.9164262 E-2 



where dist{A,B) denotes the Euchdian distance between two arbitrary points A and B 
and Q'^ represents the points of minimum distance from the station Ea to the plane V. The 
volume of the tetrahedron determined by the points Qa is proportional to the Cayley-Menger 
determinant [12] defined in terms of the lengths lij of the edges and it is given by: 



Vol 



V2 
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1 ^20 ^21 ^23 
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(37) 



which for the lengths (37) takes the value Vol = 0.7 E— 7. Therefore the metric quadruple 
(QaJab) determines a non degenerate simplex in (see for instance Saucan [IT]). This is 
equivalent to say that the metric quadruple {Qa, lab) is not congruent to any quadruple of 
points in the Euclidean plane and consequently the curvature of the wavefront surface at 
the point Qo is non vanishing. 



B. Estimation of the wavefront curvature from arrival time measurements 

Now we assume that the points {Qa}a=Q on the wavefront are directly determined through 
measurements of arrival times. Here we consider the problem of determining an approxima- 
tion of the wavefront curvature in a region far enough from the Sun (say the Earth), without 
resorting to the ray tracing method. 

An estimation of the Gaussian curvature of the wavefront surface can be obtained using 
the notion of the Wald curvature of a metric space established in the Distance Geometry (see 
[12]). that in the case of 2-dimensional manifolds agrees with the Gaussian curvature. The 
Wald curvature is determined as the limit of the embedding curvatures of metric quadruples 
isometrically embedded in surfaces of constant curvature (the Euclidean plane M."^, the 2- 
sphere or the hyperbolic space H^^). 
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FIG. 6: Schematic diagram of a wavefront plane (line P1P2) and curve (line QiQ2)- Continuous 
lines EiQi represent the light rays associated with Q1Q2 while the dashed lines EiPi represent the 
associated light rays to PiP2- 



In the hyperbolic plane of curvature — represented by the Blumenthal model 



12] . p. 19) we consider the metric quadruple {Qj,lij) defined in Subsection VA, The 



A{r) :-- 



0, 



(38) 



curvature of a hyperbolic plane on which there exists a quadruple congruent with Qj must 
fulfill both of the following conditions: 

/ 1 cosh(/oi/T) cosh(/o2/''') cosh(/o3/r)^ 

cosh(/oi/r) 1 cosh(/i2/r) cosh(/i3/r) 

cosh(/o2/'r) cosh(Zi2/r) 1 cosh(/23/''^) 

ycosh(/o3/r) cosh(/i3/r) cosh(/23/'r) 1 

and each non-zero principal minor of A{r) of order m+1 has the sign (—1)™" (see [12], p. 274). 
For small values of the curvature 1/r the determinant det (A(r)) can be approximated by a 
Taylor polynomial. Using the symbolic processor Maple and employing numerical precision 



of Digits =50 to perform the Taylor expansion we obtain that (38) may be approximated 
by the following algebraic equation for p := 1/r 



0.235 X 10"^V^° + 0.481 x 10"^^^ - 0.756 x lO'^V = 



(39) 



having only one positive real root. This approximated solution is taken as an initial guessed 



solution to solve the transcendental equation (38) by means of the command solve in 
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Maple. On the other hand, the conditions imposed on the sign of the principal minors are 
also satisfied. Therefore the Wald curvature Kw associated with the chosen quadruple {Qj} 
is 

:= = -0.16E-10. (40) 

This result gives an approximation of the total curvature of the wavefront surface under the 
assumption that locally this surface may be identified with a hyperbolic plane in which the 
quadruple considered is isometrically embedded. 

C. Scheme of the method 

In this subsection we present an outline of the construction of the method we used above 
to obtain the Wald curvature corresponding to four receiving stations located at points Ea 
and the arrival time differences Ta'- 

- Define on the wavefront the points Pa corresponding to the stations Ea- 



OQa = OEa + cr,n, (41) 

where O is the coordinate origin and n is the unit vector in the direction of the light 
rays. 

Determine the relative distances Uj between the points Qa'- 



1% = Llj + 2c(r, - Ti){OEj - OE,) ■ n + c^tj - n)'', (42) 
where Lij are the Euclidean distances between the stations Ei and Ej. 



For the point Qa so obtained, establish the nonlinear equation (38) in terms of the 
curvature of the surface. 



Solve the nonlinear equation (38) for the unknown r to obtain an estimation of the 
curvature of the wavefront in a neighborhood of the station Eq by means of the Wald 
curvature Kw of this surface modeled as a hyperbolic plane. 
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VI. CONCLUSIONS 



The ray tracing numerical method provides a useful tool for the description of spacelike 
bidimensional wavefronts within the framework of the general relativity. We have studied 
a method, based on techniques of computational geometry, that allows to estimate the 
curvature properties of the surface by making a least-squares fitting of the wavefront surface 
by a quadric surface in the neighborhood of each point of this surface. The computation of 
the light rays is carried out using an algorithm based on the Taylor method for the solution 
of differential equations and employing high arithmetic precision. Further, we have applied 
a projection at each step of the numerical integration process that allows to guarantee the 
fulfillment (at machine precision) of the isotropy condition for the tangent vector to the 
light ray. We have also studied numerically the dependence of the curvature properties of 
the wavefront surface on the value of the Eddington parameter 7. On the other hand, we 
have employed a geometric computational approach to the study of the model proposed by 
Samuel as a new general relativity test, by determining a numerical approximation of the 
volume corresponding to a tetrahedron formed by four points on the wavefront that reaches 
four receiving stations on the Earth surface. Finally, we have obtained an estimation of the 
Wald curvature for the wavefront in a vicinity of the Earth by using the differences of arrival 
time recorded at four receiving stations on the Earth. 
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